SUBROUTINE INMSOIL_NH3D_INTERFACE &
& (ix,iy,iimonth,iiday,ihour,iminu,iseco,ifcmodel,xlatit, &
& rgixiy,ratixiy, &
& SraddirUV_surf_ixiy,SraddirPAR_surf_ixiy,SraddirIR_surf_ixiy, &
& SraddifUV_surf_ixiy,SraddifPAR_surf_ixiy,SraddifIR_surf_ixiy, &
& vegixiy,albixiy,xdd2ixiy,xlakeixiy, &
& taixiy,qa,pres,precixiy,usuf,vsuf,deltaz,dt, &
& hsufixiy,z0ixiy,fcvegixiy, &
& hixiy,leixiy,cdmixiy,tsuf, &
& fmethixiy,rhoaixiy)

use MPI_VARIABLES, only : &
& nxi, nx1e, nyi, ny1e

use INMSOIL_ARRAYS

use INMSOIL_INMSOIL

use INMSOIL_PHYSPARAMS, only : &
& alb_cold_sn

use INMSOIL_NUMPARAMS, only : &
& KL, ML, MS, NT, Num_Soil, &
& Num_Veget, dznorm, dzmin, UpperLayer

implicit none

include 'soil.inc'

! Input variables
integer(4), intent(in) :: ix, iy ! Indices if gridpoint in x and y directions 
integer(4), intent(in) :: ifcmodel ! Switch for carbon models
integer(4), intent(in) :: iimonth,iiday,ihour,iminu,iseco ! time variables

real(8), intent(in) :: xlatit ! Latitude, degrees, positive North
real(8), intent(in) :: rgixiy, ratixiy ! Global shortwave radiation and longwave radiation, W/m**2
real(8), intent(in) :: SraddirUV_surf_ixiy, SraddirPAR_surf_ixiy, SraddirIR_surf_ixiy ! Direct radiation in UV, PAR and IR bands, W/m**2
real(8), intent(in) :: SraddifUV_surf_ixiy, SraddifPAR_surf_ixiy, SraddifIR_surf_ixiy ! Diffuse radiation in UV, PAR and IR bands, W/m**2
real(8), intent(in) :: vegixiy, albixiy, xdd2ixiy, xlakeixiy 
! vegixiy - fraction of gridcell covered by high vegetation
! albixiy - gridcell albedo, fraction
! xdd2ixiy - soil layer depth, m
! xlakeixiy - fraction of gridcell occupied by water bodies
real(8), intent(in) :: taixiy, qa, pres, precixiy
! taixiy - surface layer air temperature, Kelvin
! qa - specific humidity in surface air layer, kg/kg
! pres - atmospheric pressure, Pa
! precixiy - atmospheric precipitation rate, m/s
real(8), intent(in) :: usuf, vsuf ! surface air wind components in x and y directions, m/s
real(8), intent(in) :: deltaz ! height above surface at which all atmospheric variables are given, m
real(8), intent(in) :: dt ! Timestep, s
real(8), intent(in) :: hsufixiy, z0ixiy, fcvegixiy(1:nv2)
! hsufixiy - surface height above sea level, m
! z0ixiy - roughness parameter in a gridcell, m
! fcvegixiy(:) - fractions of a gridcell area occupied by different vegetation types

! Output variables

real(8), intent(inout) :: hixiy, leixiy, cdmixiy, tsuf, fmethixiy, rhoaixiy

! Local variables

real(8), parameter :: wm2_to_calm2s = 1./4.1868
real(8), parameter :: Pa_to_bars = 1.d-5
real(8), parameter :: mgCH4m2day_to_molCH4m2s = 1.d-7/(8.64*16.04)
real(8), parameter :: cm_to_m = 1.d-2, m_to_cm = 1.d+2

integer(4), parameter :: twoband = 2

real(8) :: cosz
real(8) :: fwet(1:nv2)

real(kind=8) :: solad(1:twoband)       !incoming direct solar radiation (w/m**2)
real(kind=8) :: solai(1:twoband)       !incoming diffuse solar radiation (w/m**2)
real(kind=8) :: albsod(1:twoband)      !ground albedo (direct)
real(kind=8) :: laisun(1:nv2)        !sunlite leaf area index, one-sided
real(kind=8) :: laisha(1:nv2)        !shaded leaf area index, one-sided
real(kind=8) :: parsun(1:nv2)        !average absorbed par for sunlite leaves (w/m**2)
real(kind=8) :: parsha(1:nv2)        !average absorbed par for shaded leaves (w/m**2)
real(8) :: pch
real(8), allocatable :: WSSL_work(:,:), WSSG_work(:,:)
real(8) :: rssun(1:nv2), psnsun(1:nv2), psn_unsun(1:nv2)
real(8) :: rssha(1:nv2), psnsha(1:nv2), psn_unsha(1:nv2)

real :: albsnow, albcoldsnow
real(8) :: ei, ea

integer(4) :: nsoiltype, nvegettype
integer(4) :: i, m ! Loop index

! External functions
real(8), external :: SIN_SUN, ESAT

! Assigning soil properties
nsoiltype = 0 ! nsoiltype should be a 2-d array, 
              ! but at the moment homogeneous soil
              ! over the domain is assumed (nsoiltype is a scalar)
! Assigning vegetation properties
nvegettype = 1 ! 0 means vegetation parameters are aggregated
               ! using vegetation types' coverages, should be used
               ! when fcveg is available only (e.g. when ECOCLIMAP used)

call COMSOIL( &
& KL, ML, MS, NT, NUM_SOIL, NUM_VEGET, &
& sngl(dt), DZNORM, DZMIN, sngl(xdd2ixiy), UPPERLAYER, &
& AL(1,ix,iy),DLT(1,ix,iy),DVT(1,ix,iy),ALL(1,ix,iy),DL(1,ix,iy), &
& ALV(1,ix,iy),DV(1,ix,iy),Z(1,ix,iy),T(1,ix,iy), &
& WV(1,ix,iy),WI(1,ix,iy),DENS(1,ix,iy), &
& RHOGR(1,ix,iy),WLM0(1,ix,iy),WLM7(1,ix,iy),BH(1,ix,iy),PSIMAX(1,ix,iy), &
& POR(1,ix,iy),FLWMAX(1,ix,iy),DLMAX(1,ix,iy),DZ(1,ix,iy),NSOILTYPE)

ZS(ix,iy) = hsufixiy ! surface orography, m
ST(ix,iy) = 1.
TGROLD(ix,iy) = TGRNEW(ix,iy)
QGROLD(ix,iy) = QGRNEW(ix,iy)
SNOLD(ix,iy) = SNNEW(ix,iy)
Snmelt_OLD(ix,iy) = Snmelt(ix,iy)
WSOLD(ix,iy) = WSOIL(ix,iy)
VEG(ix,iy) = vegixiy/(1. - xlakeixiy) ! the fraction of land area, covered by high vegetation
                                      ! *ixiy are fractions of the model cell
DO M = 1,ML
  T(M,ix,iy) = TSG(M,ix,iy)
  WL(M,ix,iy) = WSG(M,ix,iy)
  WV(M,ix,iy) = VSG(M,ix,iy)
  WI(M,ix,iy) = WIG(M,ix,iy)
END DO
BOLD(ix,iy) = BG(ix,iy)
AZ0P(ix,iy) = z0ixiy
PRECIP(ix,iy) = real(precixiy,4)*dt*m_to_cm ! converting meters to centimeters
!write(*,*) PRECIP(ix,iy), precixiy
!write(*,*) taixiy

! Vegetation
call VINIT( &
& sngl(taixiy),sngl(qa),sngl(pres*Pa_to_bars), & ! pressure is converted to bars
& KL, ML, MS, NT, NUM_SOIL, NUM_VEGET, &
& sngl(dt), DZNORM, DZMIN, sngl(xdd2ixiy), UPPERLAYER, &
& AZ0P(ix,iy),VEG(ix,iy),WL(1,ix,iy),ALBEDO_GR(ix,iy), &
& CCQ(ix,iy),sngl(rgixiy),TL(ix,iy),LW(ix,iy), &
& AL(1,ix,iy),DLT(1,ix,iy),DVT(1,ix,iy),ALL(1,ix,iy), &
& DL(1,ix,iy),ALV(1,ix,iy),DV(1,ix,iy),Z(1,ix,iy),T(1,ix,iy), &
& WV(1,ix,iy),WI(1,ix,iy),DENS(1,ix,iy), &
& RHOGR(1,ix,iy),WLM0(1,ix,iy),WLM7(1,ix,iy),BH(1,ix,iy),PSIMAX(1,ix,iy), &
& POR(1,ix,iy),FLWMAX(1,ix,iy),DLMAX(1,ix,iy),DZ(1,ix,iy), &
& ROOTS(1,ix,iy),WSSL(1,ix,iy), &
& WSSG(1,ix,iy),FFF(ix,iy),AMN(ix,iy),DVEG(ix,iy),RCO(ix,iy), &
& HINT(ix,iy),SVEG(ix,iy),SBARE(ix,iy),NVEGETTYPE,fcvegixiy)

! Currently snow albedo is constant
albsnow = alb_cold_sn
albcoldsnow = alb_cold_sn


!if (ix == 27 .and. iy == 27) then
!write(*,*) 'AA', &
!& sngl(usuf),sngl(vsuf), & ! wind components
!& sngl(taixiy),sngl(qa), & ! temperature and humidity
!& sngl(deltaz), &          ! atmosphere bottom level atlitude above soil
                     ! assumed for input characteristics to be
                     ! given for
!& sngl(rgixiy),  & ! shortwave radiation
!& sngl(ratixiy), & ! longwave radiation
!& VEG(ix,iy), &         ! vegetation
!& sngl(pres*Pa_to_bars),       &   ! pressure is converted to bars
!& sngl(albixiy), &    ! albedo of soil
!& albsnow, &    ! albedo of snow
!& albcoldsnow, & ! albedo of cold snow
!& T_SOIL_OUT(ix,iy), &     ! soil balanced temperature
                  ! tile fractions:
!& T_VEGET_OUT(ix,iy), &    ! tile vegetation fraction temperature
!& T_BARE_SOIL_OUT(ix,iy), & ! tile bare soil fraction temperature
!& T_WATER_OUT(ix,iy), &    ! tile water fraction temperature
!& T_SNOW_OUT(ix,iy), &     ! tile snow fraction temperature
!& Q_SOIL_OUT(ix,iy), &     ! soil humidity
!& RADIAT_OUT(ix,iy), &
!& LHFLUX_OUT(ix,iy), &     ! latent heat flux
!& SHFLUX_OUT(ix,iy), &     ! sensible heat flux
!& CDM_OUT(ix,iy), &        ! exchange coefficient
!& GSOLO_OUT(ix,iy), &      ! soil flux
!& SSBB_OUT(ix,iy), &       ! soil longwave radiation
!& PRECIP_OUT(ix,iy), &     ! precipitation accumulated at prev time step
!& METHANE_OUT(1,ix,iy), &    ! soil methane amount
!& RHOA_OUT(ix,iy), &       ! atmosphere boundary layer air density
!& EMISSIVITY_OUT(ix,iy), & ! soil emissivity     
!& KL, ML, MS, NT, NUM_SOIL, NUM_VEGET, &
!& sngl(dt), DZNORM, DZMIN, sngl(xdd2ixiy), UPPERLAYER, &
!& TSG(1,ix,iy),WSG(1,ix,iy),VSG(1,ix,iy),WIG(1,ix,iy),BG(ix,iy), &
!& WSOIL(ix,iy),WSKIN(ix,iy),WSKINNEW(ix,iy),T0(ix,iy), &
!& ROOTS(1,ix,iy),WSSL(1,ix,iy), &
!& WSSG(1,ix,iy),FFF(ix,iy),AMN(ix,iy),DVEG(ix,iy),RCO(ix,iy), &
!& HINT(ix,iy),SVEG(ix,iy),SBARE(ix,iy),ROOTSM(ix,iy), &
!& ST(ix,iy),PGR(ix,iy),TGROLD(ix,iy),QGROLD(ix,iy),RADIAT(ix,iy), &
!& PRECIP(ix,iy),WSOLD(ix,iy),SNOLD(ix,iy), &
!& ZS(ix,iy),THSOIL(ix,iy),WHSOIL(ix,iy),BOLD(ix,iy), &
!& RF1(ix,iy),RF2(ix,iy),SNMELT(ix,iy),SNMELT_OLD(ix,iy), &
!& HSNOW(ix,iy),HS(ix,iy),ES(ix,iy),TGRNEW(ix,iy),QGRNEW(ix,iy), &
!& WSNEW(ix,iy),SNNEW(ix,iy),RUNOF1(ix,iy), &
!& ELATOLD(ix,iy),HFOLD(ix,iy),PRS(ix,iy),SHORTWAVE(ix,iy), &
!& BW(ix,iy),BMELT(ix,iy),PRSOLD(ix,iy), &
!& AZ0P(ix,iy),VEG(ix,iy),WL(1,ix,iy),ALBEDO_GR(ix,iy), &
!& AL(1,ix,iy),DLT(1,ix,iy),DVT(1,ix,iy),ALL(1,ix,iy),DL(1,ix,iy), &
!& ALV(1,ix,iy),DV(1,ix,iy),Z(1,ix,iy),T(1,ix,iy), &
!& WV(1,ix,iy),WI(1,ix,iy),DENS(1,ix,iy), &
!& RHOGR(1,ix,iy),WLM0(1,ix,iy),WLM7(1,ix,iy), &
!& BH(1,ix,iy),PSIMAX(1,ix,iy), &
!& POR(1,ix,iy),FLWMAX(1,ix,iy),DLMAX(1,ix,iy),DZ(1,ix,iy), &
!& CCQ(ix,iy),SWW(ix,iy),TL(ix,iy),LW(ix,iy), &
!& EXTINCT(1,ix,iy),DENSNY(ix,iy),SN_INTRCP(ix,iy), &
!& TSRFOL(1,ix,iy),TSRF(1,ix,iy), &
!& Q_METHANE(1,ix,iy),TGRANN(1,ix,iy), &
!& FPLANT(ix,iy),FEBUL(ix,iy),FDIFF(ix,iy),FTOT(ix,iy), &
!& PLANT_SUM(ix,iy),BULL_SUM(ix,iy),OXID_SUM(ix,iy), &
!& RPROD_SUM(ix,iy),ANOX(ix,iy),WAT_TAB(ix,iy), &
!& FBL, &
!& CC(1,ix,iy), &
!& NSPIN, COND(ix,iy), PHASE(ix,iy), &
!& VKC,ANU,Z0MIN,CHC,P_M,P_T,Q_M,Q_T, &
!& ALPHA_M,ALPHA_T,BETA_M,BETA_T,GAMMA_M,GAMMA_T, &
!& ZETA_ST,A_M,A_T,B_M,B_T, &
!& WSMAX(1,ix,iy),WSL(ix,iy), &
!& HSNOLD(ix,iy), &
!& AN5(ix,iy)
!  STOP
!endif

call PBLFLX( &
& sngl(usuf),sngl(vsuf), & ! wind components
& sngl(taixiy),sngl(qa), & ! temperature and humidity
& sngl(deltaz), &          ! atmosphere bottom level atlitude above soil
                     ! assumed for input characteristics to be
                     ! given for
& sngl(rgixiy),  & ! shortwave radiation
& sngl(ratixiy), & ! longwave radiation
& VEG(ix,iy), &         ! vegetation
& sngl(pres*Pa_to_bars),       &   ! pressure is converted to bars
& sngl(albixiy), &    ! albedo of soil
& albsnow, &    ! albedo of snow
& albcoldsnow, & ! albedo of cold snow
& T_SOIL_OUT(ix,iy), &     ! soil balanced temperature
                  ! tile fractions:
& T_VEGET_OUT(ix,iy), &    ! tile vegetation fraction temperature
& T_BARE_SOIL_OUT(ix,iy), & ! tile bare soil fraction temperature
& T_WATER_OUT(ix,iy), &    ! tile water fraction temperature
& T_SNOW_OUT(ix,iy), &     ! tile snow fraction temperature
& Q_SOIL_OUT(ix,iy), &     ! soil humidity
& RADIAT_OUT(ix,iy), &
& LHFLUX_OUT(ix,iy), &     ! latent heat flux
& SHFLUX_OUT(ix,iy), &     ! sensible heat flux
& CDM_OUT(ix,iy), &        ! exchange coefficient
& GSOLO_OUT(ix,iy), &      ! soil flux
& SSBB_OUT(ix,iy), &       ! soil longwave radiation
& PRECIP_OUT(ix,iy), &     ! precipitation accumulated at prev time step
& METHANE_OUT(1,ix,iy), &    ! soil methane amount
& RHOA_OUT(ix,iy), &       ! atmosphere boundary layer air density
& EMISSIVITY_OUT(ix,iy), & ! soil emissivity     
& KL, ML, MS, NT, NUM_SOIL, NUM_VEGET, &
& sngl(dt), DZNORM, DZMIN, sngl(xdd2ixiy), UPPERLAYER, &
& TSG(1,ix,iy),WSG(1,ix,iy),VSG(1,ix,iy),WIG(1,ix,iy),BG(ix,iy), &
& WSOIL(ix,iy),WSKIN(ix,iy),WSKINNEW(ix,iy),T0(ix,iy), &
& ROOTS(1,ix,iy),WSSL(1,ix,iy), &
& WSSG(1,ix,iy),FFF(ix,iy),AMN(ix,iy),DVEG(ix,iy),RCO(ix,iy), &
& HINT(ix,iy),SVEG(ix,iy),SBARE(ix,iy),ROOTSM(ix,iy), &
& ST(ix,iy),PGR(ix,iy),TGROLD(ix,iy),QGROLD(ix,iy),RADIAT(ix,iy), &
& PRECIP(ix,iy),WSOLD(ix,iy),SNOLD(ix,iy), &
& ZS(ix,iy),THSOIL(ix,iy),WHSOIL(ix,iy),BOLD(ix,iy), &
& RF1(ix,iy),RF2(ix,iy),SNMELT(ix,iy),SNMELT_OLD(ix,iy), &
& HSNOW(ix,iy),HS(ix,iy),ES(ix,iy),TGRNEW(ix,iy),QGRNEW(ix,iy), &
& WSNEW(ix,iy),SNNEW(ix,iy),RUNOF1(ix,iy), &
& ELATOLD(ix,iy),HFOLD(ix,iy),PRS(ix,iy),SHORTWAVE(ix,iy), &
& BW(ix,iy),BMELT(ix,iy),PRSOLD(ix,iy), &
& AZ0P(ix,iy),VEG(ix,iy),WL(1,ix,iy),ALBEDO_GR(ix,iy), &
& AL(1,ix,iy),DLT(1,ix,iy),DVT(1,ix,iy),ALL(1,ix,iy),DL(1,ix,iy), &
& ALV(1,ix,iy),DV(1,ix,iy),Z(1,ix,iy),T(1,ix,iy), &
& WV(1,ix,iy),WI(1,ix,iy),DENS(1,ix,iy), &
& RHOGR(1,ix,iy),WLM0(1,ix,iy),WLM7(1,ix,iy), &
& BH(1,ix,iy),PSIMAX(1,ix,iy), &
& POR(1,ix,iy),FLWMAX(1,ix,iy),DLMAX(1,ix,iy),DZ(1,ix,iy), &
& CCQ(ix,iy),SWW(ix,iy),TL(ix,iy),LW(ix,iy), &
& EXTINCT(1,ix,iy),DENSNY(ix,iy),SN_INTRCP(ix,iy), &
& TSRFOL(1,ix,iy),TSRF(1,ix,iy), &
& Q_METHANE(1,ix,iy),TGRANN(1,ix,iy), &
& FPLANT(ix,iy),FEBUL(ix,iy),FDIFF(ix,iy),FTOT(ix,iy), &
& PLANT_SUM(ix,iy),BULL_SUM(ix,iy),OXID_SUM(ix,iy), &
& RPROD_SUM(ix,iy),ANOX(ix,iy),WAT_TAB(ix,iy), &
& FBL, &
& CC(1,ix,iy), &
& NSPIN, COND(ix,iy), PHASE(ix,iy), &
& VKC,ANU,Z0MIN,CHC,P_M,P_T,Q_M,Q_T, &
& ALPHA_M,ALPHA_T,BETA_M,BETA_T,GAMMA_M,GAMMA_T, &
& ZETA_ST,A_M,A_T,B_M,B_T, &
& WSMAX(1,ix,iy),WSL(ix,iy), &
& HSNOLD(ix,iy), &
& AN5(ix,iy),PRECIP_ACC(ix,iy))

!write(*,*) WL(ms,ix,iy), WL(ms+2,ix,iy)

DO M = 1,ML
  TSG(M,ix,iy) = T(M,ix,iy)
  WSG(M,ix,iy) = WL(M,ix,iy)
  VSG(M,ix,iy) = WV(M,ix,iy)
  WIG(M,ix,iy) = WI(M,ix,iy)
!  AT(M,ix,iy) = T(M,ix,iy)
!  AW(M,ix,iy) = WLL(M,ix,iy)
END DO
BG(ix,iy) = BOLD(ix,iy)
WSOIL(ix,iy) = WSNEW(ix,iy)
WSKIN(ix,iy) = WSKINNEW(ix,iy)
SNOW(ix,iy) = SNNEW(ix,iy)
!RUNOFF(ix,iy) = RUNOF1(ix,iy)

if (ifcmodel == 1) then
! Cosine of solar zenith angle
  cosz = SIN_SUN(iimonth,iiday,ihour,iminu,iseco,xlatit)
  fwet(1:nv2) = 0. ! Temporary solution, should be parameterized somehow
  albsod(1:twoband) = albixiy ! The same albedo for UV+PAR and IR radiation
  solad(1) = SraddirUV_surf_ixiy + SraddirPAR_surf_ixiy
  solad(2) = SraddirIR_surf_ixiy
  solai(1) = SraddifUV_surf_ixiy + SraddifPAR_surf_ixiy
  solai(2) = SraddifIR_surf_ixiy
!  write(*,*) cosz, xlatit, iimonth, iiday, ihour, iminu, iseco
  call VEGRAD(iimonth, cosz, taixiy, SNOW(ix,iy)*cm_to_m, dble(albsnow), fwet, albsod, &
  & solad, solai, laisun, laisha, parsun, parsha)

  ei = ESAT(taixiy) ! Assuming leaf temperature being equal to air temperature
  ea = ESAT(taixiy)
  pch = 1./(CDM_OUT(ix,iy)*dsqrt(usuf**2+vsuf**2)) ! resistance
  allocate (WSSL_work(1:ML,1:nv2), WSSG_work(1:ML,1:nv2))
  do i = 1, nv2
  ! Currently the same WSSG and WSSL for all vegetation types  
    WSSL_work(1:ML,i) = WSSL(1:ML,ix,iy)
    WSSG_work(1:ML,i) = WSSG(1:ML,ix,iy)
  enddo
  call STOMATA (parsun, &
  &             ei, ea, taixiy, pres, &
  &             Z(1,ix,iy), T(1,ix,iy), WL(1,ix,iy), WSSL_work, WSSG_work, &
  &             usuf, vsuf, pch, rssun, psnsun, psn_unsun)
  call STOMATA (parsha, &
  &             ei, ea, taixiy, pres, &
  &             Z(1,ix,iy), T(1,ix,iy), WL(1,ix,iy), WSSL_work, WSSG_work, &
  &             usuf, vsuf, pch, rssha, psnsha, psn_unsha)
         !write(*,*) psnsun, psnsha
  deallocate (WSSL_work, WSSG_work)
  call CBALANCE (laisun, laisha, psnsun, psnsha, dc6(1,ix,iy), &
  & fcvegixiy, dt, npp(ix,iy), c6(ix,iy), litter_in(1,ix,iy), dnpp(1,ix,iy))
  call ROMUL(litter_in(1,ix,iy), Z(1,ix,iy), T(1,ix,iy), WL(1,ix,iy), &
  & cpool_soil(1,ix,iy), npool_soil(1,ix,iy), cflux(ix,iy),dt)
endif

! Output variables
hixiy = SHFLUX_OUT(ix,iy)
leixiy = LHFLUX_OUT(ix,iy)
cdmixiy = CDM_OUT(ix,iy)*dsqrt(usuf**2+vsuf**2) ! RHOA_OUT - bug corrected
tsuf = T_SOIL_OUT(ix,iy)
fmethixiy = FTOT(ix,iy)*mgCH4m2day_to_molCH4m2s ! converting units
rhoaixiy = RHOA_OUT(ix,iy)

END SUBROUTINE INMSOIL_NH3D_INTERFACE 
